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ABSTRACT 

The origins of the hot solar corona and the supersonically expanding solar wind are still the subject of 
debate. A key obstacle in the way of producing realistic simulations of the Sun-heliosphere system is the 
lack of a physically motivated way of specifying the coronal heating rate. Recent one-dimensional models 
have been found to reproduce many observed features of the solar wind by assuming the energy comes from 
Alfven waves that are partially reflected, then dissipated by magnetohydrodynamic turbulence. However, the 
nonlocal physics of wave reflection has made it difficult to apply these processes to more sophisticated (three- 
dimensional) models. This paper presents a set of robust approximations to the solutions of the linear Alfven 
wave reflection equations. A key ingredient to the turbulent heating rate is the ratio of inward to outward 
wave power, and the approximations developed here allow this to be written explicitly in terms of local plasma 
properties at any given location. The coronal heating also depends on the frequency spectrum of Alfven waves 
in the open-field corona, which has not yet been measured directly. A model-based assumption is used here for 
the spectrum, but the results of future measurements can be incorporated easily. The resulting expression for 
the coronal heating rate is self-contained, computationally efficient, and applicable directly to global models 
of the corona and heliosphere. This paper tests and validates the approximations by comparing the results to 
exact solutions of the wave transport equations in several cases relevant to the fast and slow solar wind. 
Subject headings: interplanetary medium — MHD — solar wind — Sun: corona — turbulence — waves 



1. INTRODUCTION 

The hot, ionized outer atmosphere of the Sun is a unique 
laboratory for the study of magnetohydrodynamics (MHD) 
and plasma physics. Despite more than a half-century of 
study (Parker 1958), the basic physical processes responsi- 
ble for heating the million-degree solar corona and accelerat- 
ing the solar wind are still not known. Identification of these 
processes is important not only for understanding the origins 
and impacts of space weather (e.g., Feynman & Gabriel 2000; 
Eastwood 2008), but also for establishing a baseline of knowl- 
edge about a well-resolved star that is directly relevant to 
other astrophysical systems. 

In recent years, two general paradigms have emerged as at- 
tempts to address how both fast and slow solar wind streams 
are heated and accelerated. In the wave/turbulence-driven 
(WTD) class of models, it is generally assumed that the 
convection-driven jostling of magnetic flux tubes in the pho- 
tosphere drives wave-like fluctuations that propagate up into 
the extended corona. These waves (usually Alfven waves) 
are proposed to partially reflect back down toward the Sun, 
develop into MHD turbulence, and heat the plasma by their 
gradual dissipation (see, e.g., Hollweg 1986; Wang & Sheeley 
1991; Matthaeus et al. 1999; Suzuki & Inutsuka 2006; Cran- 
mer et al. 2007). In the reconnection/loop-opening (RLO) 
class of models, the flux tubes feeding the solar wind are as- 
sumed to be influenced by impulsive bursts of mass, momen- 
tum, and energy deposition in the low atmosphere. This en- 
ergy is usually assumed to come from magnetic reconnection 
between closed, loop-like magnetic flux systems and the open 
flux tubes that connect to the solar wind (see, e.g., Axford & 
McKenzie 1992; Fisk et al. 1999; Schwadron & McComas 
2003). 

Determining whether the WTD or RLO paradigm — or 
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some combination of the two — is the dominant cause of 
global solar wind variability is a key prerequisite to build- 
ing physically realistic models of the heliosphere. One way 
to make progress is to include either WTD or RLO processes 
in existing three-dimensional numerical simulations and com- 
pare the results with measurements. The main goal of this pa- 
per is to provide a new set of tools that allows the incorpora- 
tion of WTD physics into simulations of the Sun-heliosphere 
system. Many of the widely-applied three-dimensional mod- 
eling codes have used relatively simple empirical prescrip- 
tions for coronal heating in the energy conservation equations 
(Riley et al. 2001, 2006; Roussev et al. 2003; Toth et al. 2005; 
Usmanov & Goldstein 2006; Feng et al. 2007; Lionello et al. 
2009; Sokolov et al. 2009; Nakamizo et al. 2009; Schmit et 
al. 2009). It would also be beneficial to apply more realistic 
heating rates to more focused studies of solar wind expan- 
sion, such as two-dimensional axisymmetric models of coro- 
nal streamers (e.g., Vasquez et al. 2003; Endeve et al. 2004). 

The starting point for this work is an existing model of 
Alfven wave reflection and dissipative heating (Cranmer & 
van Ballegooijen 2005). In that model, an observationally 
constrained set of plasma parameters in a polar coronal hole 
was specified as a time-steady background state on which the 
properties of waves and turbulence were computed. Subse- 
quently, the same phenomenological wave physics was in- 
serted into a self-consistent solution of the equations of mass, 
momentum, and energy conservation (Cranmer et al. 2007). 
The only input "free parameters" to these models of coro- 
nal heating and solar wind acceleration were the photospheric 
lower boundary conditions (for the waves) and the radial de- 
pendence of the background magnetic field. For a single 
choice for the lower boundary condition, these models pro- 
duced a realistic variation of fast and slow solar wind con- 
ditions by varying only the coronal magnetic field (see also 
Cranmer 2009). There has been a great deal of other recent 
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work done to improve our understanding of Alfven wave re- 
flection and MHD turbulence as a source of coronal heating 
(e.g., Verdini et al. 2005, 2009; Verdini & Velli 2007; Rap- 
pazzo et al. 2007, 2008; Chandran & Hollweg 2009). 

Despite this progress, it has been difficult to apply the re- 
sults of these focused studies to more comprehensive models 
of the Sun-heliosphere system. Determining even the most 
basic ingredients of the theoretical coronal heating rate re- 
quires the solution of an additional set of differential equa- 
tions for the rate of Alfven wave reflection. These equations 
depend on the wave frequency, and they are inherently non- 
local in their dependence on the plasma parameters along a 
given flux tube. In other words, the amount of wave reflection 
at any given location in the corona appears to depend on an in- 
tegration over distance, and not just on the local plasma prop- 
erties. Thus, in order to compute the coronal heating from 
reflection-driven turbulence, it has been necessary to solve a 
set of differential equations for each frequency in a continuous 
power spectrum, for each flux tube of interest that fills three- 
dimensional space, and for each time step of a simulation. 

This paper presents a set of approximations that allows the 
rate of Alfven wave reflection to be computed without the 
need for computationally intensive solutions of differential 
equations. The new approximations are completely local in 
nature, in that they depend only on the plasma parameters at 
the location in the corona at which the coronal heating rate 
is to be computed. It is hoped that these approximations will 
speed up the calculation of the coronal heating rate by or- 
ders of magnitude in comparison to earlier studies. Section 
2 discusses the relevant equations and approximations. Sec- 
tion 3 compares the exact (numerically integrated) reflection 
coefficients and heating rates with those computed from the 
approximations. Section 4 describes a FORTRAN subroutine 
that has been developed to implement these approximations. 
The code is included with this paper as online-only material. 
Finally, Section 5 concludes this paper with a brief summary 
of the major results and a discussion of additional physical 
processes that can be included to improve the modeling of the 
corona and solar wind. 

2. ALFVEN WAVE REFLECTION 

This paper considers the one-dimensional variation of 
plasma parameters along a magnetic flux tube that is rooted 
in the solar photosphere and extends into interplanetary space. 
The general assumption will be that the corona and solar wind 
are in a state of steady (i.e., time independent) expansion. 
However, the heating rates discussed below may be valid un- 
der time-variable conditions as well. Throughout this sec- 
tion, the numerical examples are taken from an observation- 
ally constrained model of a polar coronal hole at solar mini- 
mum (Cranmer & van Ballegooijen 2005). 

2. 1 . The Linear Non-WKB Reflection Problem 

Alfven waves are modeled here as linear, incompressible, 
and transverse fluctuations that propagate along a magnetic 
flux tube with background field strength Bq. The wave per- 
turbations in velocity and magnetic field are denoted vj_ and 
B±. The assumption of linearity is consistent with the limiting 
case \B±/Bo\ -C 1. In general, the perturbed wave properties 
are complex quantities that vary as a function of both time t 
and heliocentric radius r. The conservation equations writ- 
ten below implicitly assume linear polarization of the waves 
along a single transverse dimension. However, this assump- 
tion does not limit the applicability of the resulting wave am- 



plitudes to other polarization states (see, e.g., Heinemann & 
Olbert 1980). 

It is convenient to express the wave amplitudes in terms of 
Elsasser (1950) variables, which are defined in velocity units 

as 

H = vx±-^=, (1) 



where p is the local mass density, z~ represents outward prop- 
agating waves, and z+ represents inward propagating waves. 
(This is a convention-dependent assignment; other papers of- 
ten use other definitions.) In a frame of reference flowing 
with the solar wind, these oscillations propagate up and down 
along the field lines with phase and group speeds equal to the 
local Alfven speed V A = B /(4irp) l l 2 . 

If the waves are propagating in only one direction along 
the field, the radial variation of their amplitude and phase can 
be described straightforwardly by defining a local wavenum- 
ber and utilizing the concept of wave action conservation 
(e.g., Jacques 1977). This limiting case is often described in 
terms of the WKB (Wentzel, Kramers, Brillouin) approxima- 
tion. However, the more general case of counterpropagating 
Alfven waves (i.e., a superposition of both Elsasser compo- 
nents) tends to require a non-WKB treatment. In this case, the 
radial evolution of the oscillation profile can no longer be ex- 
pressed with a local wavenumber. It has been known for some 
time that a spatially varying Alfven speed allows for gradual 
linear reflection (Ferraro & Plumpton 1958). This problem 
has been studied extensively in the context of solar and stellar 
winds (e.g., Hollweg 1978, 1981, 1990; Wentzel 1978; Heine- 
mann & Olbert 1980; An et al. 1990; Barkhudarov 1991; Velli 
1993; Krogulec et al. 1994; MacGregor & Charbonneau 1994; 
Orlando et al. 1996; Laitinen 2005; Verdini et al. 2005). 

For the solar models considered here, the magnitude of out- 
ward propagating waves always remains larger than the mag- 
nitude of inward (i.e., reflected) waves, and thus |z+|/|z_| < 1. 
At large heights in the corona and solar wind, it is often 
the case (for some frequencies) that the reflection is very 
inefficient, and |z+|/|z_| <C 1. It's important to note, how- 
ever, that the sharp transition region (TR) between the chro- 
mosphere and corona can act as an efficient "reflection bar- 
rier" to Alfven waves. Thus, in the photosphere and chro- 
mosphere, the reflection can be considered nearly complete 
(|z+|/|z-| ~ 1) and the fluctuations are similar in character to 
standing waves. 

The incompressible first-order conservation equations for 
mass and momentum, as well as the magnetic induction equa- 
tion, can be transformed into a pair of wave transport equa- 
tions, 



— - + (utV a )-z— 
at or 



A '\AH D 2H A J 



(2) 



where u is the solar wind speed and the signed scale heights 
are defined as Hq = p/(dp/dr) and Ha = Va/ (dV\/dr). Var- 
ious alternate ways of writing Equation (f2|i are described in 
Appendix B of Cranmer & van Ballegooijen (2005). The 
phenomenon of gradual linear reflection arises because of the 
presence of the z^/Ha term on the right-hand side. This pro- 
duces coupling between the two Elsasser variables. 

Equation (O does not contain any terms that describe the 
nonlinear damping of Alfven waves. Cranmer & van Bal- 
legooijen (2005) found that this damping does not strongly 
affect the wave amplitudes in the corona, but it may be an im- 
portant effect at larger distances in the heliosphere. Thus, it 
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appears justifiable to separate the problem of non-WKB re- 
flection from that of the nonlinear damping and coronal heat- 
ing. This is what is done in this paper. Alternately, Chan- 
dran & Hollweg (2009) presented a set of approximations of 
non-WKB reflection in which the nonlinear damping was in- 
cluded explicitly in the transport equations. It remains to be 
seen whether a combination of the approximations developed 
in this paper with those of Chandran & Hollweg (2009) will 
yield an improved description of the overall wave transport 
and dissipation. 

If the complete radial and time dependence of z+ and z~ 
were known for a given flux tube in the solar wind, it would 
be possible to compute the turbulent heating rate (see Section 
3 below). The exact solution of Equation (0, however, tradi- 
tionally requires either numerical relaxation or direct integra- 
tion up and down along the flux tube, starting at the Alfven 
critical point. This has made it difficult to incorporate an 
accurate description of reflection-driven turbulence in three- 
dimensional Sun-heliosphere simulations. 

Barkhudarov (1991) presented a dimensionless version of 
the transport equations in which the Elsasser variables are ex- 
pressed as 



z±(r,t) = G±(r)exp{i[r±(r) + wf]} 



(3) 



where the angular frequency uo (expressed in rad s" 1 ) is a real 
constant. The amplitudes G±(r) and angular phases T±{r) 
are real functions of distance along the flux tube. In order to 
determine the degree of non-WKB wave reflection, one can 
solve for two dimensionless quantities. First, a scaled ratio of 
the two amplitudes can be defined as 



* = 



u-Va \ G+ 

u+Va) gZ 



(4) 



Second, the angular phase shift between the inward and out- 
ward wave trains is defined as T = T+ - T_. Following the 
terminology of Cranmer et al. (2007), one can also define 
an effective frequency-dependent "reflection coefficient" as 
TZ= \z+\/\z~\ = G+/G-. This is the primary quantity that the 
approximations of this paper are designed to estimate. 

Barkhudarov (1991) discussed how the transport equations 
can be transformed into dimensionless conservation equations 
for the two non-WKB quantities and T. These equations are 



dV 
dr 



d^j _ (* 2 -i)cosr 
~dT ~ 2#a 

(* 2 +l)sinr 2wV A 



■n 



(5) 



(6) 



Although Barkhudarov's (1991) derivations were limited to 
the case of pure spherical expansion (i.e., Bo oc r 2 ), the above 
equations and definitions have been shown to be valid for an 
arbitrary flux-tube expansion factor (see Cranmer & van Bal- 
legooijen 2005). If the phase shift T is known, one can use 
Barkhudarov's closed-form solution to Equation (0 to solve 
for 

l-e 2w 

* = ^7 , (7) 



where 



l+e 



W(r) = / dr 



2W 



cosT 
2H7 



(8) 



In other words, if it is possible to obtain an expression for 
W in terms of radius, wave frequency, and the local plasma 



properties, then one can straightforwardly determine and 
thus K. 

The Alfven critical point (at which u = Va) is a singular 
point of the transport equations. This critical radius is denoted 
as ro. When r = ro, the wind speed and Alfven speed both have 
the identical value Vq. The ratio ^ is zero at this critical point, 
it is negative where r < ro, and it is positive where r > ro. 
Over the full range of distances, < 1. Barkhudarov (1991) 
derived the following constraint on the phase shift T at the 
singular point, 

tanr = — , (9) 



where 



du dV\ 
^ dr dr 



(10) 



It is usually the case that the Alfven speed gradient is the 
dominant contributor in the definition of \i, and one can often 
safely ignore the du/dr term above. The quantity /io is the 
value of /i at the Alfven critical point, and this (usually posi- 
tive) quantity acts as an effective "cutoff frequency" for wave 
reflection (see, e.g., Musielak et al. 1989). Waves having fre- 
quencies much lower than /l*o are strongly reflected. Waves 
having frequencies much higher than /io are reflected very 
weakly and behave similarly to WKB-like oscillations that 
obey wave action conservation. An application of L'HopitaPs 
rule gives a concise expression for the reflection coefficient at 
the critical point, which is 



Ho 



\dV A /dr\ Q 



^ 2 + Mo 



(11) 



2.2. Reflection in the Zero Frequency Limit 

In order to find approximate expressions for the amount of 
non-WKB wave reflection at heights other than the Alfven 
critical point, it is useful to examine the solutions to the di- 
mensionless wave transport equations in various liming cases. 
Numerical models show that in the limit of very low wave fre- 
quency (i.e., uj — > 0), the angular phase shift T also approaches 
zero over all radii. Thus, since cosT w 1 everywhere, Equa- 
tion ([8]) reduces to 



W = 



1 



' dlnV A 1 
dr = - In 

dr 2 



(12) 



This particularly simple solution leads to a closed-form ex- 
pression for the reflection coefficient, with 



n 7 . 



u + Va 



u-VaJ VVo + Va 



Vo-Va 



(13) 



This form of the reflection coefficient has been shown to agree 
well with the numerical solutions of Cranmer & van Balle- 
gooijen (2005) and Cranmer et al. (2007) at the lowest mod- 
eled frequencies of ~ 10~ 6 Hz. Thus, if the majority of the 
Alfven wave power in the corona and solar wind is at low 
enough frequencies (u <C /io), then Equation dT3T > provides 
the ratio of counterpropagating wave amplitudes as a function 
of only local plasma parameters (u, Va) and the velocity at the 
Alfven critical point (Vb). 

The above zero-frequency limit for the rate of linear reflec- 
tion should not be confused with a similarly named "low- 
frequency" approximation used in studies of solar wind tur- 
bulence. A series of MHD scale-separation models has been 
developed over the past few decades (e.g., Zhou & Matthaeus 
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1990; Zank et al. 1996; Matthaeus et al. 1999, 2004; Breech 
et al. 2005, 2008; Usmanov et al. 2009) in which the radial 
dependence of the power in the z+ and z~ modes is com- 
puted. The results are often expressed in terms of the nor- 
malized cross helicity, a c = (1 -1Z 2 )/(1 +1Z 2 ). In these mod- 
els the fluctuation power is assumed to be dominated by the 
lowest frequency modes. Some of these studies explicitly in- 
cluded linear wave reflection (e.g., Matthaeus et al. 1999) and 
some included only other processes such as large-scale shears, 
pickup protons in the outer heliosphere, and turbulent "dy- 
namic alignment" (for a comprehensive summary, see Breech 
et al. 2008). 

2.3. Reflection in the Infinite Frequency Limit 

Alfven waves having higher frequencies than described 
above (i.e., u> > Ho) undergo substantially weaker reflection 
than in the zero-frequency limit. In the very high frequency 
limit of u> ^> /io, the numerical models of Barkhudarov (1991) 
and Cranmer & van Ballegooijen (2005) showed that the 
phase shift T approaches an asymptotic value of -ir/2 at all 
radii. In other words, cosT — * 0, and the magnitude of W be- 
comes very close to zero everywhere, as well. This also drives 
$ and 1Z to small absolute values. At the Alfven critical point, 



cosT 



Vw 2 + // 2 



< 1 



(14) 



where the latter approximation holds for large frequencies 
(see Equation (0). It is evident from the numerical models 
of high-frequency wave reflection that Equation ( TBI serves 
reasonably well as an approximation for the entire radial de- 
pendence of the magnitude of cosT, and not merely for its 
value at ro. 

The realization that Equation ( fl4l > may be valid over all 
radii can be used to estimate W. An illustrative example is to 
examine the properties of the solar wind at large radii, where 
Va oc r~ l . This condition is usually satisfied for r > ro, and in 
some cases it is also a reasonable approximation over a few 
solar radii below ro as well. Assuming that the Alfven speed 
gradient dominates the definition of /i (i.e., Equation ([Toll), 
then pa Va/i" ~ Voro/r 2 . Applying this to Equations <[8J an d 
(fT~4T >. it becomes possible to solve for 



W 



4lo 



1 



r 2 ) 4 W 



(15) 



It was found that a slightly better approximation is to replace 
the factor of uj above with the full denominator of Equation 
(14\ . Thus, an improved approximate form for W in the limit 
of high wave frequencies is 



4y / UJ 2 + fl 2 



(16) 



where s is a dimensionless constant that usually is equal to 1, 
but sometimes needs to be set to -1 (see below). 

It should be noted that Equation (\M provides a good ap- 
proximation to the numerical solutions at large radii, but it be- 
gins to fail closer to the Sun, where the Alfven speed departs 
from an r~ l radial dependence. One problem with Equation 
( [TBI ) is that the solution for W changes sign whenever the nu- 
merator (n - no) changes sign. Most modeled radial profiles 
for VaW show at least one local maximum in the corona, and 
sometimes two or more maxima (e.g., Cranmer & van Bal- 
legooijen 2005; Evans et al. 2008). However, the numerical 
solutions for Alfven wave reflection consistently show that 
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FIG. 1 . — Approximations of the Alfven speed gradient for a polar coronal 
hole, plotted as a function of height above the solar photosphere. The mag- 
nitude of fjj is shown in regions where /i > (solid curves) and where /i < 
(dotted curves). The positive-definite estimates v (dashed curve) and Vp^/r 
(dot-dashed curve) are also shown. The filled circle denotes fiQ at the Alfven 
critical point. 

W > for r < ro, and W < for r > ro- Practically, this can 
be remedied by setting s = -1 in Equation dT6b whenever the 
condition /i < /irj occurs at radii r < ro- This is equivalent to 
taking the absolute value of Equation ( fTBI l everywhere below 
the Alfven critical point. Doing this gives values for W (and ^ 
and 1Z) that are in better agreement with the numerical mod- 
els. Nonetheless, this expression still gives rise to unphysical 
dips to W = in the narrow radial zones where /i-/io changes 
sign. The numerical models do not exhibit W = at these 
locations. 

An alternate approximation is to replace the derivatives in 
the definition of /i with a positive-definite expression that in- 
volves only the plasma parameters themselves. Experimen- 
tation with a range of functional forms led to the definition 
of 

„= ^ (17) 

(r+R Q )(r-R Q ) 

where R Q is the solar radius, and it can be seen that v « \x for 
large radii r 3> Rq- This expression can be substituted for /i 
in Equation (fT6] i to form an alternate definition for W, 



V—Vo 
Wocv = 

Wuj 2 + v 2 



(18) 



where vq is the value of v at the Alfven critical point. Figure 
1 shows the radial dependence of /i and v, as well as the even 
simpler approximation \i « VaA used above, for the coronal 
hole model of Cranmer & van Ballegooijen (2005). Note that 
v « | | nearly everywhere, but v remains continuous and pos- 
itive at the locations in the corona where [i changes sign. The 
simpler expression Va/V disagrees significantly with the mag- 
nitudes of both fi and v at low heights in the corona. 

In practice, an even better high-frequency estimate of 1Z 
was found by taking the arithmetic average of the resulting re- 
flection coefficients that come from Equations ( TToT ) and ( fl8l ). 
As can be seen below, there is a moderate amount of "spik- 
iness" in the numerically computed reflection coefficients at 
the radii where fj, changes sign. Using only Equation ( [ToT l 
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would have overestimated that effect, and using only Equa- 
tion < fT~8T > would have underestimated it. 

2.4. Bridging for Intermediate Frequencies 

The two extreme limits of low frequencies (Section 2.2) and 
high frequencies (Section 2.3) can be combined to produce a 
more complete estimate for the coupled radius and frequency 
dependence of the reflection coefficient. The most robust way 
of doing this is to first compute the separate values of 1Z for 
each of the above limiting forms of W, and then combine them 
in the following way: 



n = nl: 



(19) 



The dimensionless bridging exponent e needs to be a function 
of both radius and wave frequency. When < e < 1 , the above 
functional form gives a resulting value of 1Z that is intermedi- 
ate between the low and high frequency limiting cases. Some 
experimentation found the optimal definition 



1 



e = 



l + (w /w) 2 



(20) 



where ujq = (v + vq)/2. At heights where ui is much smaller 
than the local value of loq, the exponent e w and Equation 
( fT9l is dominated by the low-frequency limit lZ zero . On the 
other hand, when uj too, the exponent e approaches 1 and 
the bridging relation is dominated by the average of the two 
alternate forms of the high-frequency limit for 1Z. 

What is the behavior of the "bridging frequency" loq! At 
small radii (i.e., in the corona) where the local value of v 
is large compared to uq, the bridging frequency loq is domi- 
nated by that local value. However, at large radii above the 
critical point, where v <C vq, the bridging frequency is not al- 
lowed to decrease below Vo/2. This behavior forms a bound- 
ary between the high and low frequency regions (in radius- 
frequency space) that matches what is seen in the numerical 
solutions. 

It should be noted that the above estimates are intended to 
apply to the corona and the solar wind, and not to the pho- 
tosphere and chromosphere. For the latter (low-temperature, 
high-density) regions that sit below the sharp TR, a relatively 
safe approximation would be to simply assign R» 1. In 
any case, it is likely that in the chromosphere, other sources 
of heating — such as the entropy gain at shocks formed by 
the steepening of acoustic waves — are more important than 
Alfven waves (e.g., Cranmer et al. 2007). 

2.5. Estimating Properties of the Alfven Critical Point 

In one-dimensional models of the plasma conditions along 
a specified magnetic flux tube, the location of the Alfven criti- 
cal point can be found rather easily. In that case, the numerical 
values of Vo, /xq, and vq would also be known. However, in 
multi-dimensional MHD simulations, in which the conserva- 
tion equations are solved either by discretization or by spec- 
tral methods, it may not be feasible to calculate these quanti- 
ties for each point in space and time. Doing so would require 
computationally intensive integrations up and down along in- 
dividual magnetic field lines, at each time step. Since the main 
purpose of this paper is to eliminate the need for similar kinds 
integrations of the non-WKB equations, it would be advanta- 
geous to find approximate ways of computing the properties 
of the Alfven critical point — even when all one knows are the 
properties at a location far from this point. Thus, in order to be 



able to apply the techniques developed in Sections 2.2-2.4 to 
as wide a range of models as possible, this subsection presents 
an "optional" method to estimate ro and Vo when these are not 
known a priori. 

For a steady-state solar wind, the condition of mass flux 
conservation demands that the quantity pu/Bo remain con- 
stant along any flux tube. Substituting in the definition of 
the Alfven speed, this condition implies that the density at the 
Alfven critical point is determined uniquely to be 



(21) 



po = constant = p \ — 

Va 



where all quantities on the right-hand side are evaluated at 
any arbitrary location along the flux tube. It is interesting 
that a value for po can be computed even if neither ro nor 
Vq are known for that flux tube. The ratio p/po is useful as a 
definitive probe of whether the current location in the corona 
or solar wind is below (p/po > 1) or above {p/po < 1) the 
Alfven critical point. 

In the numerical models of Cranmer & van Ballegooijen 
(2005) and Cranmer et al. (2007), the value of Vb is always 
intermediate between the instantaneous values of u and Va- 1 
Thus, a parameterization of the form 



(22) 



was found to be useful, where < a < 1. Some trial-and- 
error experimentation led to a density-dependent fit for the 
exponent, with 



1 



l+0.3(p/po) 



0.25 



(23) 



Figure 2 shows the radial dependence of the estimated value 
of Vb for the coronal hole model of Cranmer & van Ballegooi- 
jen (2005). This approximation gives a roughly constant mag- 
nitude for Vo throughout most of the corona and solar wind. 
However, there is likely to be room for improvement in the 
choices for the numerical constants in Equation ( 1231 . This 
parameterization should be tested and refined by comparing 
with additional models of the corona and solar wind. Once 
Vo is known, a reasonable approximation for the radius of the 
Alfven critical point is 



ro 



rV A 
Vo 



(24) 



This is an exact expression over the range of heights at which 
Va oc r~ l , but it appears to also provide a decent order-of- 
magnitude estimate for ro at other heights as well. 

Although ro and Vo can be estimated using the methods out- 
lined in this section, the derivatives required for computing po 
are likely to be less amenable to robust approximation. Thus, 
when ro and Vo are estimated in this way, we suggest that 
Equation ( fT6l ) be avoided for the high-frequency limit of the 
reflection coefficient. In this case, the bridging law given by 
Equation ( [T9b should be replaced by a simpler version, 



T? = 7? 1_£ 7? c 



(25) 



1 This does not apply to the chromospheric and photospheric regions below 
the TR. At those low heights, both the local values of u and Va often dip below 
Vb. In these regions, however, the methods outlined in Sections 2.2-2.5 do 
not need to be used and TZ = 1 is not a bad approximation. 
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FIG. 2. — Radial dependence of th e estimated speed Vq at the Alfven critical 
point computed from Equation (22) (solid curve). This is compared with the 
radial dependence of Va (dashed curve) and u (dotted curve) for this coronal 
hole model. The estimate for Vq is shown only above the sharp chromosphere- 
corona transition region at (r/R@) — 1 ~ 0.003. 

3. RESULTS 

This section contains a detailed comparison between the 
numerically computed non-WKB wave properties (i.e., exact 
solutions to Equations © and ©) and the results of the ap- 
proximations described in Section 2. In Sections 3.1 and 3.2, 
the background coronal properties shown are those of Cran- 
mer & van Ballegooijen (2005). In Section 3.3, the approxi- 
mate heating rates are compared to those computed by Cran- 
mer et al. (2007) for models of a range of source regions of 
fast and slow solar wind streams. 

3.1. Reflection Coefficients 

The non-WKB reflection equations were solved numeri- 
cally using an adaptive-stepsize version of the fourth order 
Runge-Kutta algorithm (e.g., Press et al. 1992). The base- 
line model consists of a grid of 350 frequency points and 
5457 radial points. The frequency points are evenly spaced in 
logw. The minimum and maximum wave periods are 0.001 
and 1000 hours, which correspond to maximum and minimum 
frequencies of 2.77 x 10 _1 and 2.77 x 10~ 7 Hz respectively. 
The radial grid extends from the mid-chromospheric "merg- 
ing height" of strong-field flux tubes (r sa 1. 00086 Rq) to a 
heliocentric distance of 4 AU (r f=s 860/?©), and the radially 
varying grid separation is described by Cranmer & van Balle- 
gooijen (2005). 

Figure 3(a) shows the result of numerically integrating the 
transport equations to compute TZ as a function of radius for 
each "monochromatic" frequency in the grid. As described in 
the caption, the contours in all three panels denote constant 
values of TZ between 3 x 10~ 5 and 0.9. Figure 3(b) displays 
the result of using the approximations of Sections 2.2-2.4 to 
compute TZ from Equation dl9) . For Figure 3(b), the known 
values of r = 9.698/? Q , V = 660.1 km s _1 , and hq = 1.184 x 
10~ 4 s" 1 were applied to the approximation equations. Figure 
3(c), however, shows the result of using the estimation method 
of Section 2.5 to compute ro and Vq, and to approximate the 
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FIG. 3. — Contour plot of the reflection coefficient TZ as a function of 
height and wave frequency. Exact numerical results (a) are compared with 
results of applying approximations of Sections 2.2-2.4 (b), as well as the 
more approximate results of estimating critical point conditions using Section 
2.5 (c). The 3 panels use identical definitions for the contour levels. From 
light to dark (i.e., from bottom to top along the right-hand side), the solid 
contours represent values of 71 of 0.9, 0.3, 0.1, 0.03, 0.01, 0.003, 0.001, 
3 X 10" 4 , 10" 4 , and 3 X 10~ 5 . Dotted lines in (b) and (c) show the bridging 
frequency (u>o(r)/2Tr) for the two approximations. 

value of TZ from Equation d25l l. 

Overall, the similarities between the three panels of Figure 
3 appear to outweigh the differences. At large heights (r > 
5Rq) the high-frequency limiting approximations for TZ pro- 
duce excellent agreement with the numerical solutions. The 
full radial dependence at very low frequencies (w/27r < 10~ 6 
Hz) is well reproduced by the results of Section 2.2. There 
are some discrepancies between the three panels at heights in 
the low corona (between 0.003 and 0.1 Rq) at low and in- 
termediate frequencies. However, these discrepancies involve 
mainly the positions of the contours labeling TZ = 0.1 and 0.3, 
and the approximations here do not show more than a fac- 
tor of two difference from the numerical solutions. The ef- 
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fects of n changing sign are apparent at heights of about 0.02 
and 0.5 Rq for the highest frequencies, but these "spikes" are 
smoothed over for frequencies lower than about 0.01 Hz. The 
average of the two high-frequency approximations (Equation 
([191 )) produces a small cusp at these heights, which is an ade- 
quate representation of the mean radial behavior of 7Z at these 
frequencies. 

In order to apply the computed reflection coefficients to 
a model of turbulent dissipation and heating, one needs to 
understand the properties of the continuous power spectrum 
Pa(oj) of Alfvenic fluctuation energy in the corona and solar 
wind. In other words, one needs to know how much is con- 
tributed by each frequency (i.e., each "row" in Figure 3) to 
the total turbulent energy. Although there are observational 
constraints on the frequency dependence of P\(u) at the so- 
lar photosphere (Cranmer & van Ballegooijen 2005) and there 
are direct in situ measurements at distances greater than 60 Rq 
(Tu & Marsch 1995; Goldstein et al. 1997), there are no mea- 
surements of the Alfvenic frequency spectrum in the extended 
corona and inner solar wind. Unfortunately, the heights with- 
out measured power spectra appear to be the most important 
for the specification of the rates of heating and momentum 
deposition in the models. 

In this paper, the frequency-weighting of the modeled re- 
flection coefficients will be presented for two empirically 
based choices for the shape of Pa(w). Figure 4 displays both 
spectra, each of which has been normalized such that 



duj P a (oj) = 1 



(26) 



For simplicity, the normalized spectrum is assumed to retain 
its shape as a function of height. This assumption has served 
reasonably well in existing models of solar wind acceleration 
(e.g., Cranmer et al. 2007), but in reality there must be some 
kind of spectral evolution with increasing distance from the 
Sun. A model of this evolution can be straightforwardly com- 
bined with the approximations of Section 2, but that is beyond 
the scope of this paper (see, however, Tu & Marsch 1995; 
Laitinen 2005; Verdini et al. 2009). 

The solid curve in Figure 4 shows the power spectrum of 
total (kinetic plus magnetic) wave energy from the model 
of Cranmer & van Ballegooijen (2005) taken at the base of 
the corona. This model was constrained at the photospheric 
lower boundary by a measured power spectrum of the ki- 
netic motions of G-band bright points. These bright points 
represent thin magnetic flux tubes that undergo both random 
walks, in response to convective granulation, and rapid hori- 
zontal "jumps" that appear to be the result of sporadic merg- 
ing and fragmenting. The spectrum at the coronal base was 
computed as the result of a non-WKB model of kink-mode 
and Alfven-mode transport in the chromosphere and transi- 
tion region. The dashed curve in Figure 4 shows a power-law 
spectrum reminiscent of Kolmogorov (1941) hydrodynamic 
turbulence. The lower and upper limits in frequency were 
obtained from the non-WKB reflection models of Verdini & 
Velli (2007), who studied the implications of a oj~ 5 / 3 power 
law spectrum on the coronal heating rate. 

For comparison, Figure 4 also shows a measured power 
spectrum of Alfvenic fluctuations in the low corona (Tom- 
czyk & Mcintosh 2009) from a series of Doppler images 
made by the Coronal Multi-channel Polarimeter (CoMP) at 
the Sacramento Peak Observatory. It should be noted that the 
coronal regions that dominated the measured spectrum were 
closed loops in a coronal active region. Because a turbulent 
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FIG. 4. — Normalized Alfven wave power spectra for the model of Cranmer 
& van Ballegooijen (2005) (solid curve) and for a Kolmogorov (1941) power 
law with frequency limits from Verdini & Velli (2007) (dashed curve). Also 
shown, with arbitrary normalization, is the measured coronal power spectrum 
from Tomczyk & Mcintosh (2009) (dotted curve). 



cascade behaves somewhat differently in open and closed re- 
gions (e.g., Rappazzo et al. 2008; Cranmer 2009), there may 
not be a good reason to assume that this spectrum would ex- 
ist in the source regions of the solar wind. There also may 
be frequency-dependent line-of-sight integration effects that 
change the shape of the spectrum. It is nonetheless interest- 
ing that this measured spectrum appears to combine the over- 
all power-law behavior of the Kolmogorov (1941) curve with 
a hint of the high-frequency convective resonance features in 
the Cranmer & van Ballegooijen (2005) spectrum. 
The spectrum-weighted reflection coefficient (1Z) is defined 



{1l)\r) 



f du;P A (uj)TZ 2 (uj,r) 
J dujP A (uj) 



(27) 



where the square of 1Z is used because the power spectrum 
is an energy density quantity and 1Z is a ratio of amplitudes. 
Equation d27l i was used in the solar wind acceleration mod- 
els of Cranmer et al. (2007). A slightly different technique 
was applied by Cranmer & van Ballegooijen (2005), who per- 
formed the weighting on the energy densities of the z+ and 
Z- fluctuations separately, and then computed their ratio after- 
wards. The resulting heating rates from the two techniques 
are not significantly different from one another. 

In Figure 5, the radial dependences of the weighted reflec- 
tion coefficients (1Z) are shown for the exact and approximate 
calculations. Note that 1Z (for all frequencies) is set to 1 at 
heights below the TR. This is obviously not an exact repre- 
sentation of the numerically integrated value, but — as can be 
seen below in Figure 6 — the resulting heating rate is insensi- 
tive to these differences. 

2 If the number of oscillations along the line of sight at any one time scales 
with wavelength, then higher frequencies would undergo more line-of-sight 
"Doppler cancellation" than lower frequencies. Thus, the intrinsic coronal 
spectrum would have to be flatter than the measured spectrum. This could 
imply a closer resemblance to the Cranmer & van Ballegooijen (2005) spec- 
tram than is apparent in Figure 4. 
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FIG. 5. — Weighted reflection coefficients, shown as a function of radial 
distance, for (a) the Cranmer & van Ballegooijen (2005) spectrum and (b) 
the truncated Kolmogorov (1941) spectrum. In each panel, exact numeri- 
cal results (solid curves) are compared with approximate results computed 
with known rg and Vq (dashed curves) and the approximate results computed 
with local estimates for ro and Vo from Section 2.5 (dotted curves). Mea- 
sured ranges for fast-wind reflection coefficients from Helios and Ulysses are 
shown as individual error bars. 

Because the Kolmogorov spectrum is dominated by the 
lowest frequencies, the reflection is dominated by the zero- 
frequency limit of Section 2.2. Thus, the comparison in Fig- 
ure 5(b) between the exact result and the approximation com- 
puted with known values of rg and Vo shows nearly exact 
agreement. The curves in Figure 5(a) are dominated by higher 
frequencies and show more of a relative discrepancy between 
the exact and approximate expressions. For both panels, the 
reflection coefficients (TZ) computed with the estimated val- 
ues of ro and Vo (i.e., the dotted curves) show a significantly 
more "smoothed out" radial dependence than the other two 
sets of curves. This gives rise to smoother radial variations in 
the heating rates. Although this smoothing represents a de- 
parture from the exact non-WKB results, it is probably not a 
major problem. Complete models of coronal heating and solar 
wind acceleration must also contain heat conduction, and no 
matter what fine radial structure may exist in the heating rate 
itself, the existence of conduction gives rise to a similar kind 
of smearing of the thermal energy along the magnetic field. 

Figure 5 also shows Helios and Ulysses measurements in 
the fast solar wind that have been processed from the Elsasser 



energy densities presented by Bavassano et al. (2000). The 
upper and lower limits of the error bars reflect the spread in 
individual data points, each of which represented one-hour 
datasets. The frequency range of the fluctuations sampled by 
these measurements spanned only about an order of magni- 
tude, from 3 x 10~ 4 Hz to 4 x 10~ 3 Hz. As can be seen in 
Figure 4, this range of frequencies is where the two theo- 
retical spectra have comparable power to one another. It is 
the frequencies outside this narrow range that give rise to the 
significant differences between the modeled sets of curves in 
Figures 5(a) and 5(b). Contributions from higher frequencies 
drive down the reflection coefficient in Figure 5(a), and con- 
tributions from lower frequencies drive it up in Figure 5(b). 
Thus, it is not surprising that the measurements of power at 
intermediate frequencies give reflection coefficients that fall 
between the two sets of curves. 

3.2. Turbulent Hea ting Rates 

The adopted phenomenological rate of coronal heating is an 
expression for the total energy flux that cascades from large to 
small scales. It is constrained by the properties of the fluctua- 
tions at the largest scales, and it does not specify the exact ki- 
netic means of dissipation once the energy reaches the small- 
est scales. Dimensionally, this is similar to the rate of cas- 
cading energy flux derived by von Karman & Howarth (1938) 
for isotropic hydrodynamic turbulence. The full form, which 
takes into account various MHD effects, is 



2turb — P^turb 



ztz + +zlz. 



(28) 



(Hossain et al. 1995; Zhou & Matthaeus 1990; Matthaeus et 
al. 1999; Dmitruk et al. 2001, 2002; Breech et al. 2008). The 
individual components of Equation (|2B1 are described below. 

The absolute values of the spectrum- weighted Elsasser vari- 
ables are denoted Z+ and Z_. Strictly speaking, in a model of 
non-WKB wave reflection the power in both components can 
be specified only by combining the derived reflection coeffi- 
cient TZ with a specification of the absolute power spectrum 
of the fluctuation energy. However, the assumption that the 
total wave power varies in accord with straightforward wave 
action conservation has been shown to be reasonable, even in 
interplanetary space where TZ is not small (Zank et al. 1996; 
Cranmer & van Ballegooijen 2005). The examples shown be- 
low utilize this assumption in order to compute Z+ and Z_. 

Under wave action conservation, the product of a (modi- 
fied) energy flux and the transverse area of the flux tube re- 
mains constant as a function of radial distance. For disper- 
sionless Alfven waves, this is equivalent to maintaining a con- 
stant energy flux per unit magnetic field strength, or 



F 
So" 



(u + V A ) 2 U A 
VaB () 



= constant , 



(29) 



where the wave energy density Ua = pv\ (see, e.g., Jacques 
1977). Measurements of v_l can be made by analyzing the 
nonthermal Doppler broadening of coronal emission lines. At 
r 1.1 Rq, measurements made by the SUMER instrument 
on SOHO corresponded to a range of perpendicular velocity 
amplitudes vj_ w 50-60 km s" 1 (Banerjee et al. 1998; Landi & 
Cranmer 2009). Using the plasma parameters from the coro- 
nal hole model discussed above, these amplitudes are consis- 
tent with a range of values for F /Bq between about 5 x 10 4 
and 8 x 10 4 erg cm" 2 s" 1 G _1 . At r 1.5/?©, lower and upper 
limit measurements made by the UVCS instrument on SOHO 
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gave a range of amplitudes v_l ~ 100-140 km s 1 (Esser et 
al. 1999). These are consistent with a range of F/Bq values 
similar to those at the lower height: 6 x 10 4 to 10 5 erg cm -2 
s" 1 G" 1 . 

The empirically constrained models of Cranmer & van Bal- 
legooijen (2005) were used as another way to put limits on 
the most likely wave amplitudes in the corona. These models 
included wave dissipation due to the presence of turbulence, 
so the "constant" F/Bq actually decreases slightly with in- 
creasing distance. In those models, the numerical values of 
F/Bq ranged between about 2 x 10 4 and 9 x 10 4 erg cm" 2 s" 1 
G -1 . For the remainder of this paper, an intermediate con- 
stant value of 5 x 10 4 erg cm -2 s" 1 G -1 is assumed. With 
that, the full radial dependence of Ua can be computed, and 
the spectrum-averaged Elsasser amplitudes can be determined 
from 



4t/ A 



Z + = (TZ)Z- 



(30) 



for use in Equation I 

The perpendicular length scale L± is an effective transverse 
correlation length of the turbulence for the largest eddies. The 
models presented here use a standard assumption that the cor- 
relation length scales with the transverse width of the mag- 

-1/2 

netic flux tube; i.e., that L±_ oc B Q (see also Hollweg 1986). 
Cranmer & van Ballegooijen (2005) found that the following 
normalization produced a reasonable combination of heating 
and wave dissipation: 



L ± (r) « 1 L55 Mm , 
VMr) 



(31) 



where Bo is measured in Gauss. However, the self-consistent 
models of Cranmer et al. (2007) worked best by reducing the 
normalization by about a factor of four; i.e., by replacing the 
factor of 1 1.55 in the above by the smaller value 2.876. With 
no convincing evidence to choose between these two values, 
the example models shown below use a compromise value of 
6.0 Mm for this normalization. 

The quantity £ tl „b in Equation ( |28l l is an efficiency factor 
that attempts to account for regions where the turbulent cas- 
cade may not have time to develop before the fluctuations are 
carried away by the wind. Cranmer et al. (2007) estimated 
this efficiency factor to scale as 



1 



turb 



l+feddyAref)" 



(32) 



where the two timescales above are f e ddy, a nonlinear eddy cas- 
cade time, and f re f , a timescale for large-scale Alfven wave re- 
flection. The value n = 1 was chosen for the exponent based on 
a range of analytic and numerical turbulence models (Pouquet 
et al. 1976; Dobrowolny et al. 1980; Matthaeus & Zhou 1989; 
Dmitruk & Matthaeus 2003; Oughton et al. 2006). Equa- 
tion ( l32t quenches the turbulent heating when f e dd y ^> ? re f> i- e -> 
when the Alfven waves want to propagate away much faster 
than the cascade can proceed at a given location. Cranmer et 
al. (2007) defined the reflection time as ? re f = Vl^ ' ^a|- For 
the purposes of this paper, the simpler approximation f re f = v~ l 
was found to give an equivalent result (see Equation (fTTIl). 
This latter definition is used in the results presented below. 
The eddy cascade time is given by 



^eddy 



Lj_ V37T 
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FIG. 6. — Turbulent heating rate per unit mass (Qt U rb/p) for the coronal hole 
models. The curves shown in panels (a) and (b), and the line plotting styles 
in each panel, are the same as in Figure 5. The gray regions show empirically 
constrained heating rates (see the text). 

where the Alfven Mach number Ma = u /Va and the numeri- 
cal factor of 3ir comes from the normalization of an assumed 
shape of the turbulence spectrum (see Appendix C of Cranmer 
& van Ballegooijen 2005). 

Figure 6 shows the computed heating rates for the various 
spectrum-weighted cases discussed above. Rather than plot 
<2turb itself, Figure 6 shows the heating rate per unit mass 
Qtuib/p in order to more clearly indicate which heights re- 
ceive the most heating on a particle-by-particle basis. As in 
Figure 5, the differences between the exact and approximate 
non-WKB results are relatively minor. 

Each panel of Figure 6 also shows observational constraints 
on the heating rate in the fast solar wind. The gray region in 
the corona (i.e., at heights between 0.2 and 6 Rq above the 
surface) delineates the lower and upper bounds on a set of pa- 
rameterized heating rates taken from a range of papers (Wang 
1994; Hansteen & Leer 1995; Allen et al. 1998; Cranmer et 
al. 2007). These are heating rates that are required for these 
one-dimensional coronal models to be able to produce realis- 
tic fast solar wind conditions. Figure 2(b) of Cranmer (2004) 
illustrates some of the individual heating functions that were 
collected into this set. The gray region in interplanetary space 
(i.e., at heights greater than ^60 Rq) illustrates the range of 
total (Qp + Qe) heating rates determined empirically by Cran- 
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mer et al. (2009) from Helios and Ulysses plasma measure- 
ments in the fast wind. 

Note that the heating rates computed from the high- 
frequency-dominated spectrum of Cranmer & van Ballegooi- 
jen (2005) — shown in Figure 6(a) — tend to be in better agree- 
ment with the empirical constraints than do the heating rates 
computed from the low-frequency-dominated spectrum, in 
Figure 6(b). This is consistent with the results of Cranmer 
et al. (2007), in which a similar high-frequency-dominated 
spectrum was used to produce reasonably successful models 
of both fast and slow solar wind streams. Note also that both 
sets of model results in Figure 6 overestimate the measured 
heating in interplanetary space. This is likely to be the result 
of modeling the radial dependence of Z± with wave action 
conservation. In the models of Cranmer & van Ballegooijen 
(2005), the turbulent dissipation produced a factor of 4 reduc- 
tion in v_l at 1 AU, from 144 km s" 1 (undamped) to 35.2 km 
s" 1 (damped). This would give about a factor of 4 3 ' 2 = 8 de- 
crease in the heating rate far from the Sun. Inserting this into 
the present models would push the curves in Figure 6(a) to the 
bottom of the empirical range, and it would push the curves 
in Figure 6(b) down to the top of the empirical range. 

3.3. Fast and Slow Wind Source Regions 

Although the numerical examples shown above were com- 
puted for the fast solar wind that emerges from a polar coro- 
nal hole, the approximations developed in this paper are not 
meant to be limited to only that type of solar wind. Figure 7 
shows heating rates computed not only for coronal holes, but 
also source regions of slow wind at solar minimum (i.e., the 
legs of quiescent equatorial streamers) and at solar maximum 
(i.e., active regions). These models are self-consistent numer- 
ical calculations of the non-WKB reflection, turbulent heat- 
ing, and wind acceleration that were described by Cranmer et 
al. (2007). The quiescent streamer model corresponds to the 
"last" open field line that originates at a colatitude of 29.7° 
in the two-dimensional axisymmetric magnetic field model of 
Banaszkiewicz et al. (1998). The active region model has an 
added component to its magnetic field strength that was pa- 
rameterized by Cranmer et al. (2007) as having an exponen- 
tial height dependence of BAe~ (r ~ R °^ h , where = 50 G and 
/i = 0.07/? Q . 

Figure 7(a) shows the radial dependence of the wind speed 
u for these three models. Note that the active region model has 
an outflow speed of order 100 km s" 1 in the low corona. This 
appears to correspond with recent measurements of active re- 
gion outflows of similar magnitude made by the Extreme- 
ultraviolet Imaging Spectrometer (EIS) on Hinode (see, e.g., 
Harra et al. 2008). Figure 7(b) presents the heating rates per 
unit mass (Qt m b/p) that were computed self-consistently in 
the Cranmer et al. (2007) models. Figure 7(c) shows the heat- 
ing rates as computed under the approximations of this paper. 

For the active region slow wind model, two curves are 
shown in Figure 7(c). The lower curve was computed us- 
ing the standard form for the dimensionless turbulent effi- 
ciency factor £ tur b given by Equation (l32l . The upper curve 
was computed by not using this efficiency factor (i.e., by as- 
suming £ tul b = 1). Clearly, the strong peak in the heating rate 
that is seen in Figure 7(b) is reproduced much more satis- 
factorily when the efficiency factor is not used. The other 
two models — for the coronal hole and streamer cases — did 
not exhibit as strong a relative difference in the heating rate 
when £ tur b was set to unity. It should be noted that this fac- 
tor was not used in the earliest published versions of Equation 
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FIG . 7 . — Comparisons of plasma parameters for source regions of fast and 
slow solar wind: i.e., polar coronal holes (solid curves), quiescent equatorial 
streamers (dashed curves), and active region field lines connected to the solar 
wind (dotted curves). Numerical results for (a) outflow velocity, and (b) heat- 
ing rate per unit mass from Cranmer et al. (2007) are compared with heating 
rates computed using the HEATCVB code (c). 



d28l i (e.g., Hossain et al. 1995), nor is it used in other more 
recent studies (e.g., Chandran & Hollweg 2009). Thus, this 
term should still be viewed as an approximate attempt to take 
account of "wave escape" effects and not the final word in 
that story. Future improvements in modeling the turbulence- 
driven solar wind should involve finding a more robust way 
of expressing how turbulence is quenched in regions where 
waves escape more rapidly than they can cascade and dissi- 
pate. 

4. NUMERICAL CODE 

A brief FORTRAN 77 subroutine called HEATCVB has 
been developed to implement the approximations and other 
expressions given in Sections 2 and 3. The source code is 
included with this paper as online-only material, and it is 
also provided, with updates as needed, on the author's web 
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page. 3 There are only four required inputs to this subroutine: 
the radial distance r, the wind speed u, the mass density p, 
and the magnetic field strength |Z?o|- The primary output of 
the subroutine is the heating rate Q tur b at the location defined 
by the input parameters. All quantities are assumed to be in 
cgs/Gaussian units. 

In addition to the four required input parameters, the 
HEATCVB subroutine has three optional input parameters: 
the radius of the Alfven critical point ro, the speed at this lo- 
cation (i.e., Vb = u(ro) = V\(ro)), and the local Alfven wave 
velocity amplitude vj_. If the input values of ro or Vo are less 
than or equal to zero, the code will recompute them using the 
estimation method described in Section 2.5. If the input value 
of vj_ is less than or equal to zero, it will be recomputed using 
wave action conservation (i.e., Equation (|29|i). This approx- 
imation appears to be reasonably valid in the corona (where 
most of the heating occurs), but it overestimates the Alfven 
wave amplitudes in interplanetary space. 

Another derived quantity that is computed and output by 
HEATCVB is the local turbulent dissipation rate of the waves, 
7 = fiturb/(2t/A)- This quantity is specified in units of s" 1 and 
can be used when solving for departures from wave action 
conservation. The factor of 2 in the above definition assumes 
the standard convention that 7 is the damping rate for the wave 
amplitude. If instead the damping rate for wave energy is re- 
quired, 7 = Qturb/t^A- Cranmer & van Ballegooijen (2005) 
showed how this kind of wave dissipation is necessary for 
agreement between measurements of Alfven wave amplitudes 
in the corona and in the heliosphere (see also Roberts 1989; 
Verdini & Velli 2007). 

For simplicity, the HEATCVB code does not utilize the p 
gradients described in Section 2.3. Thus, the bridging be- 
tween the low-frequency and high-frequency limits for TZ 
is computed using Equation d25l ) instead of Equation $1% . 
The weighting of TZ over the Alfvenic power spectrum is 
performed using the high-frequency-dominated spectrum of 
Cranmer & van Ballegooijen (2005). As described above, 
the use of this particular choice of the spectrum worked well 
for the self-consistent solar wind models of Cranmer et al. 
(2007). The integration over the normalized power spectrum 
is discretized into bins separated by factors of two in fre- 
quency space (i.e., "octaves"). The wave power in each bin 
i is summed into weighting factors fi, such that Equation d27| > 
is approximated by 

(K) 2 « £>ft 2 (u>/)- (34) 

i 

To make the above compatible with Equation d27l i. the 
weights fi are defined such that they sum to 1. The 
HEATCVB code uses 17 bins that span more than four or- 
ders of magnitude in frequency space. The code also contains 
several checks to ensure that the input parameters are within 
reasonable bounds for the solar atmosphere and solar wind. 
All major steps in the algorithm are described in comments 
within the source code. 

The HEATCVB routine contains several tests to evaluate 
whether a full calculation of the non-WKB reflection is war- 
ranted. Whenever u < 0, the code avoids this calculation and 
sets (TZ) = 1 . The general assumption in this case is that the 
local plasma is in a closed-field region. The condition u = 
could correspond to a hydrostatic coronal streamer, and u < 
could signal the existence of transient downflows similar to 

3 http://www.cfa. harvard.edu/~scranmer/ 



those observed with visible-light coronagraphs (Wang et al. 
1999). In flux tubes where both footpoints are anchored to the 
solar surface, Alfven waves are believed to propagate up and 
down with nearly equal intensities in the two Elsasser com- 
ponents (see, e.g., Rappazzo et al. 2008). Also, whenever 
p>2x 10~ 13 g cm -3 the local plasma conditions are judged 
to be "chromospheric." In this case, the code sets (TZ) = 1 be- 
cause it assumes the location to be modeled is below the sharp 
reflection barrier of the TR. This density criterion may also 
be triggered for coronal mass ejections (CMEs). These re- 
gions may be similar to the above case of coronal loops, since 
the field lines are closed in CME plasmoids, and the shocked 
regions in front of the plasmoids are generally believed to 
have both footpoints rooted to the solar surface (e.g., Lin et 
al. 2003). However, in CMEs the assumption of time-steady 
wave action conservation may not give an accurate value for 
the local fluctuation amplitude v±. 

5. CONCLUSIONS 

The primary aim of this paper has been to develop and test 
a set of approximations to the solutions of the equations of 
non-WKB Alfven wave reflection. These solutions are de- 
signed to be applied to modeling the plasma heating in open- 
field regions of the solar corona and the solar wind. Two in- 
dependent approximate expressions for the reflection coeffi- 
cients TZ were developed in the limiting cases of extremely 
low wave frequencies (Section 2.2) and extremely high wave 
frequencies (Section 2.3). These, together with a technique 
to estimate the radius and speed at the Alfven critical point, 
provide a robust approximation for the necessary ingredients 
of the phenomenological turbulent heating rate. The resulting 
radial dependence of TZ and the heating rate Q m ±, were shown 
to agree well with exact solutions of the non-WKB transport 
equations in several cases relevant to the fast and slow solar 
wind. 

Because the approximations presented above do not de- 
pend on computationally intensive integrations along flux 
tubes, they make it much easier to insert more realis- 
tic (wave/turbulence heating) physics into three-dimensional 
models of the corona and heliosphere. An important first step 
is to perform "testbed" simulations for specific time periods 
in which large amounts of empirical data are available. Such 
times include the Whole Sun Month in 1996 (Galvin & Kohl 
1999) and the Whole Heliosphere Interval in 2008 (Gibson et 
al. 2009). These testbed models can help determine whether, 
for example, the WTD or RLO paradigms for solar wind ac- 
celeration are better representations of the actual physics (see 
Section 1). These models will also be key to assessing and 
validating real-time predictive models of the plasma condi- 
tions in the heliosphere. 

A significant remaining unknown quantity in the above 
modeling methodology is the shape of the frequency spec- 
trum of Alfven waves in the corona. The HEATCVB code 
utilizes an empirically constrained spectrum from the work 
of Cranmer & van Ballegooijen (2005). This particular form 
of the spectrum was shown by Cranmer et al. (2007) to pro- 
duce reasonably successful self-consistent models of both fast 
and slow solar wind streams, so it appears to be a good first 
approximation. However, this can be replaced easily with 
other forms of the power spectrum if better constraints be- 
come available. It is also likely that the shape of the power 
spectrum depends on radial distance and cannot be specified 
universally for all values of r (e.g., Verdini et al. 2009). Once 
such a radial dependence is specified, however, the approxi- 
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mations developed in this paper can still be used to estimate 
the turbulent heating at any location along the open field lines. 

In order to make further progress, it will be important to in- 
clude other physical processes in WTD-type models of coro- 
nal heating and solar wind acceleration. In addition to non- 
WKB reflection, there are other ways that the energy in out- 
ward propagating waves can be tapped to produce inward 
propagating waves. Large-scale shear motions in the helio- 
sphere have been suggested as a source of kinetic energy that 
could go into reducing the overall cross helicity of the turbu- 
lence (e.g., Roberts et al. 1992; Zank et al. 1996; Matthaeus et 
al. 2004; Breech et al. 2008; Usmanov et al. 2009). Kinetic in- 
stabilities in collisionless regions of the corona and solar wind 
can give rise to the local generation of high-frequency in- 
ward propagating waves (Isenberg 2001 ; Isenberg et al. 2009). 
Also, outward propagating Alfven waves of sufficient am- 
plitude may become unstable to nonlinear processes — such 
as "parametric decay" — which give rise to enhanced inward 
propagating wave activity (e.g., Goldstein 1978; Jayanti & 
Hollweg 1993; Lau & Siregar 1996; Ofman & Davila 1998; 
Suzuki & Inutsuka 2006; Hollweg & Isenberg 2007). 

Although the phenomenological heating rate given in Equa- 
tion (l28l was motivated by the results of many analytic and 
numerical studies, there is still some disagreement about its 
robustness and applicability. Chandran et al. (2009) found 
that the exact dependence of the dissipation rate on Z + and 
Z_ may depend on the specific generation mechanism(s) of 
inward propagating waves. Also, the evolution of the cor- 
relation length L± with radial distance should be coupled to 
the overall evolution of Z± and not simply scaled with the 
flux-tube width as in Equation d3~TT l (see, e.g., Matthaeus et al. 
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Finally, it will be important for future models to take ac- 
count of the multi-fluid nature of coronal heating and solar 
wind acceleration. A large-scale description of the energy flux 
injected into the turbulent cascade is needed in order to model 
its eventual kinetic dissipation (and the subsequent preferen- 
tial energization of electrons, protons, and heavy ions). Even 
in the strongly collisional "low corona," there are likely to 
be macroscopic dynamical consequences of the partitioning 
of energy between protons and electrons (see, e.g., Hansteen 
& Leer 1995; Endeve et al. 2004). Non-WKB Alfven wave 
reflection also affects the energy and momentum coupling be- 
tween protons and other ions in the solar wind (Li & Li 2007, 
2008). Including these effects can lead to concrete predictions 
for measurements to be made by space missions such as So- 
lar Probe (McComas et al. 2007) and Solar Orbiter (Marsden 
& Fleck 2007), as well as next-generation ultraviolet corona- 
graph spectroscopy that could follow up on the successes of 
the UVCS instrument on SOHO (Kohl et al. 2006). 
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